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Abstract 

A  method  for  automatically  aligning  consecutive  datasets  of  large,  two-dimensional  multi-tile 
electron  backscatter  diffraction  (EBSD)  scans  with  high  accuracy  was  developed.  The  method 
involves  first  locating  grain  and  phase  boundaries  within  search  regions  containing  overlapping 
data  in  adjacent  scan  tiles,  and  subsequently  using  cross-correlation  algorithms  to  detennine  the 
relative  position  of  the  individual  scan  tiles  which  maximizes  the  fraction  of  overlapping 
boundaries.  Savitzky-Golay  filtering  in  two  dimensions  was  used  to  estimate  the  background, 
which  was  then  subtracted  from  the  cross-correlation  to  enhance  the  peak  signal  in  samples  with 
a  high  density  of  interfaces.  The  technique  was  demonstrated  on  datasets  with  a  range  of  grain 
boundary  densities.  The  equations  were  implemented  as  enhancements  to  a  recently  published 
open  source  code  for  automated  stitching  of  multi-tile  datasets. 

I.  Introduction 

Electron  backscatter  diffraction  is  an  advanced  scanning  electron  microscope  (SEM)  based 
technique  for  characterizing  the  local  microstructure  and  crystallography  of  materials.  The 
technique  was  commercialized  in  the  early  1990’s  and  has  seen  widespread  adoption  in  a  variety 
of  fields,  most  notably  geology  and  materials  science.  Advances  in  hardware  [Schwarzer,  2008, 
Schwarzer,  2010]  and  software  [ED AX,  Inc.,  2007;  Bruker  Nano  GmbH,  2010]  technology  are 
pushing  acquisition  speeds  in  excess  of  1  kHz  with  further  improvements  on  the  horizon. 
Naturally,  this  leads  to  the  collection  of  larger  datasets,  especially  when  combined  with  fast, 
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precise  stage  movements  [Shiveley,  2011].  This  creates  two  new  issues  that  must  be  overcome: 
(1)  how  to  best  combine  the  individual  datasets,  or  tiles,  and  (2)  how  to  best  analyze  these 
datasets  that  contain  tens  of  millions  of  data  points  with  file  sizes  that  can  exceed  1  GB  in  raw 
ASCII  format.  The  latter  issue  is  being  tackled  by  commercial  vendors  who  are  taking  advantage 
of  64-bit  architecture  in  their  latest  software  versions,  which  allows  access  to  much  larger  blocks 
of  memory.  The  first  issue,  however,  has  received  considerably  less  attention.  In  previous  work, 
an  open  source  code  [Pilchak,  2011]  was  developed  to  facilitate  combination,  or  stitching,  of 
two-dimensional,  multi-tile  datasets  in  a  variety  of  fonnats  allowing  for  precise,  user-defined 
positioning  of  the  tiles.  The  software  is  especially  well  suited  for  datasets  with  coarser  step  sizes 
mapping  large  areas  where  precise  grain  morphology  is  not  an  absolutely  critical  attribute  of  the 
analysis.  The  code  also  produces  generally  acceptable  results  on  datasets  with  fine  (sub¬ 
micrometer)  step  sizes  where  grain-level  morphology  is  of  interest,  but  the  inaccuracies  inherent 
to  screw-driven  SEM  stages  with  rotary  encoders  inevitably  lead  to  subtle  variations  in  the  true 
tile  to  tile  overlap  which  introduces  artifacts  into  the  seams  between  adjacent  tiles.  This  is  most 
frequently  observed  when  the  stage  has  to  travel  a  long  distance  to  return  to  the  beginning  of  a 
new  row.  In  this  work,  we  expand  on  our  prior  software  and  demonstrate  procedures  for  high 
accuracy,  automated  combination  of  EBSD  data  from  adjacent  tiles  in  a  multi-tile  dataset  using 
the  information  contained  within  the  individual  tiles. 

II.  Background 

Cross-correlation  and  convolution  are  widely  used  for  signal  and  image  processing  in  many 
scientific  disciplines.  Because  of  their  similarities,  there  is  often  confusion  regarding  the  two 
processes  and  so  they  are  introduced  with  a  simple,  yet  relevant  example.  Cross-correlation  and 
convolution  are  both  linear  and  shift-invariant  operations,  meaning  that  the  same  operation  is 
applied  at  every  point  in  the  discrete  representation  of  the  image  (or  signal).  Both  involve 
performing  an  entry-wise  multiplication  of  two-dimensional  arrays  which  result  in  the  creation  of 
a  new  array.  Consider  the  rectangular  image  /  and  the  filter  G,  shown  in  Figure  1(a)  and  (b), 
respectively.  These  are  binary  representations  of  the  right  most  and  left  most  50  pm  and  the 
entire  vertical  size  (188  pm,  at  0.5  pm  resolution)  of  two  adjacent  tiles  in  a  multi-tile  dataset  of  a 
hot  pack  rolled  and  recrystallized  nickel-base  superalloy.  The  cross-correlation  of  these  two 
arrays  is  defined  as: 
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N  N 

G°I(x,y)=^j  ’y^,G/llI(y+k)(v+n 

l=-N  k=-N  (1) 

where  k  and  /  represent  the  row  and  column  indices  of  pixels  in  the  filter  G,  (x,y)  are  the 
coordinates  of  the  pixels  in  /  and  2N+1  is  the  number  of  elements  in  G.  The  results  of  the  cross¬ 
correlation  of  Figure  1(a)  and  Figure  1(b)  is  shown  in  Figure  1(c).  Mathematically,  convolution 
differs  only  in  that  the  filter  is  mirrored  both  horizontally  and  vertically,  i.e. 

G*I{x,y)=fjfjGklI(x_kXy_l) 

I=-Nk=-N  (2) 

The  mirrored  filter  is  shown  in  Figure  1(d),  and  the  resulting  convolution  of  Figure  1(a)  and  (b), 
or  identically,  the  cross-correlation  of  Figure  1(a)  and  (d),  is  shown  in  Figure  1(e).  The  insets  in 
Figure  1(c)  and  (e)  show  the  location  around  the  highest  intensity  at  a  magnification  of  300 
times.  It  is  obvious  that  the  cross-correlation  produces  a  strong  peak  intensity  above  the 
background  noise.  The  convolution,  on  the  other  hand,  produces  mainly  background  noise  with 
numerous  local  maxima  around  the  strongest  intensity.  This  is,  of  course,  because  /  (the  filter)  is 
mirrored  and  therefore  can  never  be  perfectly  matched  with  G,  the  image,  unless  it  is  symmetric. 
Thus,  the  signal  in  the  convolution  is  primarily  related  to  the  density  of  grain  boundary  points 
and  the  unifonnity  of  their  distribution  whereas  cross-correlation  is  a  direct  process.  The 
background  noise  is  essentially  due  to  the  fact  that  there  are  a  nonuniform  number  of  pixels  that 
contribute  to  the  solution  to  the  cross-correlation  and  convolution  for  a  given  (x,y)  coordinate. 
Thus,  spurious  grain  boundary  noise  is  maximized  at  the  central  position  where  the  search 
regions  are  fully  overlapped  and  tends  to  decrease  toward  the  periphery  of  the  solution.  In  fact,  it 
was  noted  that  the  shape  of  the  background  was  similar  to  that  of  the  instantaneous  area  of  the 
overlapping  regions  of  G  and  /  when  each  pixel  in  the  cross-correlation  or  convolution  were 
calculated  as  evidenced  by  Figure  1(f). 

This  simple  example  demonstrates  why  cross-correlation  is  used  for  pattern  recognition,  template 
matching  and  image  registration.  Practical  applications  include  automatic  alignment  of  SEM 
stages  during  three-dimensional  serial  sectioning  experiments  [Groeber,  2006;  Uchic,  2006]  and 
for  measurement  of  surface  strains  [Kammers,  2001;  Walley,  2012],  a  technique  commonly 
referred  to  as  digital  image  correlation.  Convolution,  on  the  other  hand,  is  commonly  used  to 
enhance  specific  features  of  an  image.  Uses  include  reducing  noise  prior  to  segmentation  and 
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feature  extraction  [Simmons,  2009].  In  the  standard  method  of  indexing  EBSD  patterns,  a 
‘butterfly’  convolution  filter  is  used  to  enhance  the  peaks  in  Hough  space  associated  with 
specific  diffracting  (hkl)  planes  [Schwarzer,  1997].  In  image  processing,  edge  detection  filters 
are  used  to  enhance  grain  or  phase  boundaries  prior  to  segmentation  [Simmons,  2009].  In  this 
work,  we  demonstrate  the  use  of  standard  cross-correlation  algorithms  to  automatically  combine 
multi-tile  EBSD  datasets  collected  with  a  combination  of  beam  scans  and  stage  movements. 

III.  Procedure 

First,  we  briefly  describe  the  overall  methodology  for  automated  stitching  followed  by  a  more 
detailed  treatment  of  the  individual  steps.  In  this  work,  previously  published  source  code 
[Pilchak,  2011]  was  modified  to  automatically  determine  the  linear  shifts  necessary  to  bring  each 
tile  of  a  multi-tile  EBSD  dataset  into  registration.  In  this  formulation,  the  data  points  in  each  tile 
must  be  collected  on  a  square  grid,  though  the  tiles  themselves  need  not  have  the  same 
dimensions  in  the  X  and  Y  directions.  The  major  structure  of  the  code,  that  is,  the  file  naming 
convention,  the  sequence  in  which  the  tiles  are  opened  and  the  rules  for  deciding  the  location  of  a 
particular  tile  (i.e.  first  tile,  remaining  tiles  in  the  first  row,  each  first  tile  of  a  new  row  of  data, 
and  internal  tiles)  was  unchanged.  Additional  functions  were  added  as  necessary  to  perform 
automated  stitching  as  described  below. 

Each  tile  was  opened  sequentially  and  treated  as  the  current  tile.  The  neighboring  tile(s)  relevant 
for  determining  its  placement  were  also  opened.  The  data  points  in  each  tile  that  belonged  to  the 
so-called  search  region  were  extracted  and  the  grain  and/or  phase  boundaries  within  it  were 
identified  resulting  in  binary  arrays  similar  to  those  shown  in  Figure  1(a)  and  (b).  A  schematic  of 
the  various  cases  encountered  during  stitching  are  shown  in  Figure  2.  Depending  on  which  case 
is  encountered,  the  size  and/or  orientation  of  the  search  region  may  change.  For  border  tiles,  such 
as  those  in  the  first  row  or  column,  tile  placement  depends  on  only  one  neighbor.  On  the  other 
hand,  tile  placement  depends  on  two  neighbors  for  internal  tiles.  In  this  case,  the  search  regions 
from  the  two  neighboring  tiles  are  stitched  using  the  previously  calculated  shifts  to  create  an 
inverted  L  shaped  search  region. 
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Next,  cross-correlation  was  performed  on  the  grain  boundary  arrays  which  revealed  the  position 
that  maximized  the  number  of  overlapping  grain/phase  boundaries.  In  datasets  with  a  high 
density  of  boundaries,  a  two-dimensional  Savitzky-Golay  filter  was  applied  to  the  result  of  the 
cross-correlation  to  provide  an  estimate  of  the  background.  This  result  was  subtracted  from  the 
cross-correlation  to  increase  the  signal  to  noise  ratio  in  order  to  identify  the  linear  shifts 
necessary  to  align  the  tiles.  These,  in  turn,  were  used  to  update  the  (x,y)  coordinates  for  each 
pixel  in  the  current  tile  before  appending  it  to  the  output  file.  The  approach  described  above  was 
implemented  in  MATLAB®  and  used  the  OdfPf  function  set  distributed  by  the  Deformation 
Processes  Laboratory  at  Cornell  University  [Dawson,  201 1]  to  perfonn  misorientation 
calculations  quickly. 

A.  Identification  of  Grain  and  Phase  Boundaries 

Grain  and  phase  boundaries  were  identified  using  a  [2x2]  kernel  in  a  fashion  similar  to  standard 
edge  detection  in  image  processing  except,  instead  of  comparing  grayscale  intensity,  the 
misorientation  angle  between  each  pixel  and  its  neighbors  in  the  kernel 

"0  f 

-1  ^  (3) 

were  used  to  detect  the  grain  boundaries  or  “edges”.  The  kernel  was  implemented  as  a  sliding 
window  operation  which  involved  calculating  the  misorientation  angle  between  the  pixel  beneath 
the  0  and  its  neighbors  identified  by  1  ’s  in  the  kernel.  If  the  scalar  misorientation  angle  in  any  of 
the  neighboring  bins  exceeded  a  user-defined  threshold,  typically  in  the  range  5°  to  10°,  the  pixel 
beneath  the  0  was  assigned  a  value  of  1  thereby  identifying  that  it  was  adjacent  to  a  grain 
boundary.  Similarly,  the  integer  identifying  the  phase  associated  with  each  point  in  the  kernel 
was  also  compared  with  its  neighbors  also  assigned  a  1  if  the  phases  were  different.  The  result  of 
these  operations  was  a  binary  array  in  which  the  pixels  adjacent  to  grain  boundaries  are  equal  to 
1  and  grain  interiors  are  equal  to  0.  The  method  for  locating  grain  boundaries  is  demonstrated  in 
Figure  3. 

Calculations  of  misorientation  can  be  performed  efficiently  using  the  quaternion  representation 
of  rotations.  Unit  quaternions  are  4-index  vectors  which  have  the  fonn  Q  =  [q0,q]  =  [cos(co/2),  n 

sin(w/2)]  where  qO  and  q  are  the  scalar  and  vector  components,  respectively,  and  satisfy  the 
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relationship  Q=  \ql  +  q]  +  q\  +  ryj]12  ■  Omega  (0)  represents  the  rotation  angle  while  the  unit 

vector  n  represents  the  rotation  axis  of  the  crystallite  with  respect  to  the  sample  reference  frame 
(we  adopt  Nye’s  convention  [Nye,  1985]).  The  quaternion  misorientation  in  the  fundamental 
zone,  which  yields  the  minimum  misorientation  angle  between  two  quaternions,  Qa  and  Qb,  can 
be  determined  by  finding  the  minimum  rotation  angle  corresponding  to  the  product  of  the  Q} 
and  Q“ : 

A Q  =  QaQb  =  [do, a  qo,Ba  +  qA'qBH,  qo,AqBa  -  qo,Ba*qA  -  Qa  *  qBfl]  (4) 

where  the  superscript  a  indicates  the  use  of  all  symmetrically  equivalent  orientations  of  Qb. 

B.  Cross-correlation  of  search  regions 

Once  the  grain  boundaries  in  the  search  region  have  been  identified,  the  next  step  is  to  detennine 
the  relative  position  of  the  arrays  that  produces  the  maximum  amount  of  boundary  overlap.  The 
basic  fonn  of  the  calculation  is  similar  for  the  three  types  of  tiles,  namely,  those  in  the  first  row 
of  data,  the  first  column  of  data  and  internal  tiles  having  neighbors  both  above  and  to  the  left. 

The  code  was  originally  developed  with  the  cross-correlation  process  being  performed  in  the 
spatial  domain  before  ultimately  converting  to  frequency  domain  calculations  for  significantly 
reduced  computational  costs. 

While  the  search  region  arrays  containing  the  grain  boundaries  are  discrete  functions,  it  is  still 
possible  to  take  advantage  of  the  convolution  theorem  [Bracewell,  1999]  which  states  that  the 
Fourier  transfonn,  F,  of  a  convolution  (or  correlation)  of  two  functions  is  equivalent  to  the  entry- 
wise  product  of  the  Fourier  transforms  of  the  individual  functions: 

F{LoC}=F  {L}  •  (F  {C}  )*  (5) 

where  L  and  C  are  two-dimensional  binary  arrays  representing  the  grain  boundaries  in  adjacent 
search  regions  in  the  left  and  current  tiles  in  the  first  row  of  data  and  the  asterisk  denotes  the 
complex  conjugate.  Thus,  cross-correlation  according  to  Eq.  (1),  where  L  was  fixed  and  C  is  the 
filter,  can  be  performed  quickly  to  determine  the  appropriate  shifts  that  align  L  and  C  which 
occurs  at  the  location  of  the  peak  intensity  in  the  solution  to  Eq.  (5).  The  offset  in  pixels  (x0,y0) 
can  be  recovered  from  the  cross-correlation  with: 

Xo  /max  ~  N C'  L  (6) 
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(7) 


y0  kmax  -  Nrx 

where  (kmax,lmay)  is  the  location  of  the  peak  intensity  in  F{ L  o  C}  and  Nr, l  and  /Vc,l  are  the 
number  of  rows  and  columns  in  the  search  region  of  the  left  tile.  These  translate  to  shifts  of 

A'C  —  XL  +  Xsize  -  X(yTstep  (8) 

VC  —  V0Tstep  (9) 

in  the  spatial  domain  where  (xc,  vc)  are  the  updated  coordinates  of  the  current  tile.  It  is  evident 
that  the  new  x-coordinates  depend  on  those  of  the  prior  tile,  L,  but  the  v-coordinates  do  not  for 
the  first  row  of  tiles. 

Analogous  equations  were  formulated  for  the  case  when  the  current  tile  (C)  is  the  first  tile  in  a 
new  row: 

F{UoQ=F  {U}  ■  (F  {C})*  (10) 

and  the  offset  is  recovered  and  applied  with: 

X0  /peak  "  Nc,  U  (1 1) 

Vo  kpeak  ”  F r  L  (12) 

A(  A\,.Vsk.p  (13) 

Vc  —  Vu  F  Tsize  F  V„  f step  (14) 

Internal  tiles,  with  neighbors  above  and  to  the  left,  were  treated  in  a  different  manner.  The  left 
and  upper  tiles  are  first  opened  and  their  boundaries  registered  in  a  new  array,  N,  using  the 
values  previously  calculated  to  create  an  inverted  L  shape  with  the  remaining  part  of  the  array 
being  padded  with  zeros.  The  origin  of  the  left  tile  search  region  in  N  is  located  at: 


i  =  K,  u- 


yend, v  -  (y  end,  L  ^ 'size ) 


while  the  origin  of  the  search  region  for  the  upper  tile  is  at: 


k  =  NcX  - 


^end,  L  ^Fend,  U  Fslz(, ) 


where  (xe„d,  yend)  represent  the  coordinate  (in  micrometers)  of  the  lower  right  coordinate  of  either 
the  left  (L)  or  upper  (U)  tiles.  Then,  the  boundaries  in  a  similar  inverted  L  shaped  search  region 
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of  the  current  internal  tile  were  extracted.  These  two  arrays  were  correlated  in  the  frequency 
domain, 

F{N°C}  =  F{n}-(F{C})* 

and  the  offset  was  recovered, 


*0  ^peak  ”  NC)C  (20) 

yo  ^peak  “  Nr>c  (21) 

and  applied  in  the  spatial  domain  as  a  function  of  position  with  respect  to  the  left  tile: 

•^C  Xs[Zq  ~  X o^-step  (22) 

yc  ~~  yL  +  ^size  +  y0Ystep  (23) 


C.  Background  subtraction 

Background  noise  caused  by  improperly  overlapped  grain  boundaries  was  evident  throughout  the 
correlations.  Generally,  the  background  was  sufficiently  small  that  it  did  not  obscure  the  peak  in 
the  cross-correlation,  but  samples  with  high  densities  of  boundaries  required  the  background  to 
be  estimated  and  subtracted  from  the  correlation  before  the  correct  peak  was  able  to  be 
automatically  identified.  In  early  approaches,  each  point  in  the  correlation  was  divided  by  the 
instantaneous  area  of  the  overlapping  regions  of  the  search  regions  at  the  time  that  pixel  was 
calculated.  While  appropriate  for  some  datasets,  area  nonnalizing  caused  creation  of  some 
artificial  peaks  near  the  edges  of  the  correlation,  where  the  divisor  was  small,  and  on  others  with 
more  dense  boundaries.  In  addition,  the  shape  of  the  background  was  dependent  on  the 
morphology  of  the  microstructure,  and  so  a  more  robust,  adaptive  method  was  pursued  that 
would  estimate  the  background  directly  from  each  solution  to  the  correlation  functions. 

To  this  end,  Savitzky  and  Golay’s  [Savitzky,  1964]  smoothing  algorithm  was  used  to  estimate 
the  background  of  the  two-dimensional  cross-correlation  function.  This  filter  is  known  for 
preserving  local  maxima  and  minima  that  otherwise  may  be  obscured  by  simple  Gaussian,  mean 
and  median  smoothing  filters.  The  smoothed  signal  was  subtracted  from  the  real  signal  to 
produce  a  background  corrected  cross-correlation.  To  obtain  the  smoothed  signal,  an  odd  number 
of  consecutive,  equidistant  data  points  (up  to  25  in  the  original  formulation,  though  the  equations 
have  been  generalized  [Madden,  1978])  were  fit  to  a  polynomial  of  the  form: 
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(24) 


p(x )  -  w0  +  wtx  +  w2x2  + ...  +  wnx"  =  ^  WmXm 

m= 0 

where  w  were  smoothing  coefficients  using  a  least  squares  fit.  Eq.  (24)  is  typically  applied  as  a 
sliding  window  operation  in  order  to  smooth  an  entire  dataset.  Thus,  each  pixel  becomes  the 
center  element  of  the  sliding  window  and  the  value  in  the  smoothed  dataset  is  equal  to  the 
solution  to  Eq.  (24)  evaluated  at  the  same  position  once  the  coefficients  are  known.  Savitzky  and 
Golay  noted  that  rather  than  perfonning  brute  force  least  squares  regression,  the  coefficients  can 
be  calculated  easily  using  standard  matrix  methods, 

Aw  =  c  (25) 

where  w  is  a  column  vector  of  the  coefficients,  c  is  a  column  vector  of  the  values  neighboring 
pixels  in  the  kernel  (centered  on  pixel  kf)  and  A  is  the  design  matrix  Au  =  kl.  By  multiplying  both 
sides  by  the  transpose  of  A, 

AtAw  =  Atc  (26) 

it  can  be  seen  that 

w  =  (AtA)'1Atc  (27) 

These  equations  are  solved  twice,  first  row-wise,  then  column-wise  on  the  already  row-smoothed 
data,  in  order  to  produce  the  two-dimensional  background  necessary  to  subtract  from  the 
correlations,  Eqs.  (5),  (10)  and  (19),  if  necessary. 

IV.  Results 

In  this  section,  we  plot  solutions  to  several  of  the  equations  presented  above  for  three  datasets  of 
varying  complexity.  The  first  is  a  simple,  single  phase  Ti-7A1  (wt%)  alloy  with  equiaxed  grains 
of  the  order  of  the  tile  size.  The  second  is  a  2  row  by  2  column  subset  of  tiles  from  a 
polycrystalline  nickel-base  superalloy  prepared  via  hot  pack  rolling  and  recrystallization  heat 
treatment.  The  third  is  an  a  +  P  forged  Ti-6A1-4V  alloy  with  large  microtextured  regions.  In  the 
Ti  alloy  literature,  microtextured  regions  [Sinha,  2007],  also  called  macrozones  [Gennain,  2007], 
are  known  as  clusters  of  equiaxed  a  grains  that  have  similar  c-axis  orientation  in  the  hep  lattice 
but  otherwise  appear  recrystallized  /  equiaxed  in  the  morphological  sense.  These  are  formed  by 
the  geometric  breakdown  of  the  similarly  oriented  platelets  in  a  colonies  [Bieler,  2002]  that  were 
formed  during  a  prior  P  annealing  step  in  the  thermomechanical  processing  sequence.  These 
datasets  were  chosen  because  they  have  a  range  of  phase  and  grain  boundary  densities. 
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Following  each  of  the  detailed  analyses,  the  code  is  demonstrated  on  the  entire  data  sets  from  the 
Ti-6A1-4V  and  nickel  superalloy  samples. 


Table  1.  EBSD  data  acquisition  parameters  for  the  three  demonstration  datasets. 


Material 

Tile  size,  pm 

SEM  mag. 

Number  of  tiles 

(rows,  columns) 

Step  size, 

pm 

Search  range, 

pm 

Ti-7A1 

115  x  115 

800x 

30  (5,  6) 

0.2 

15 

Ti-6A1-4V 

270  x  270 

350x 

650  (25,  26) 

2 

20 

Ni-15A1- 

5Cr+C,B,Zr 

187  x  187 

500x 

44  (4,  11) 

0.5 

40 

A.  Coarse  grain,  single  phase  a-titanium 

The  coarse  grained  titanium  alloy  serves  as  an  instructive  example  because,  as  evident  in  Figure 
4(a),  only  one  triple  point  and  one  grain  boundary  appear  in  the  search  regions  for  the  left  and 
current  tiles.  Consequently,  there  is  essentially  no  background  noise  and  the  contributions  of 
individual  grain  boundaries  are  evident  as  lines  of  increased  intensity  in  the  solution  to  Eq.  (5), 
Figure  5(c).  These  lines  have  similar  slope  and  orientation  to  the  grain  boundary  that  causes 
them.  As  expected,  the  maximum  value  of  the  cross-correlation  is  located  at  the  intersection  of 
these  lines  indicating  that  this  is  the  location  of  best  possible  match.  There  was  very  little 
background  noise  in  this  coarse  grained  material,  so  the  maximum  in  the  cross-correlation  was 
sufficient  to  cover  the  offset.  For  thoroughness,  however,  the  Savitzky-Golay  estimated 
background  is  shown  in  Figure  5(d)  while  Figure  5(e)  shows  the  difference  between  the  cross¬ 
correlation  and  the  smoothed  background.  In  the  cross-correlation,  the  maximum  value  is 
approximately  185  while  the  adjacent  regions  are  ~0  (dark  blue  regions)  and  ~60  (light  blue 
regions  corresponding  to  grain  boundaries)  giving  a  signal  to  noise  ratio  of  approximately  3:1.  In 
contrast  the  background  corrected  peak  intensity  is  ~48  while  the  background  noise  is  ~5,  a  ratio 
of  nearly  10: 1.  The  results  of  this  cross-correlation  are  shown  in  Figure  6.  The  grain  boundaries 
corresponding  to  only  the  overlapping  parts  of  left  and  current  tiles  when  the  cross-correlation  is 
maximized  are  shown  in  Figure  6(a)  and  Figure  6(b),  respectively,  while  the  entrywise  product 
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of  these  arrays  appears  in  Figure  6(c).  The  inverse  pole  figure  map  in  Figure  6(d)  reveals  good 
matching  with  no  disconnects  or  seams  created  on  the  grain  boundaries  in  the  stitched  dataset. 

B.  Polycrystalline  nickel-base  superalloy 

A  2  row  by  2  column  dataset  was  combined  in  this  example  to  show  the  case  for  stitching  an 
internal  tile.  The  cross-correlation  for  the  first  two  tiles  appeared  as  the  first  example,  Figure 
1(c).  The  presence  of  a  slight  background,  whose  shape  was  similar  to  the  area  of  overlapping 
regions  at  each  point  in  the  cross-correlation  was  noted,  but  the  peak  in  the  correlation  was  still 
more  intense  and  therefore  easily  located.  The  overlapping  parts  of  the  left  and  current  tile  arrays 
at  the  point  where  the  correlation  is  maximized,  and  the  entrywise  (element-by-element)  product 
of  these,  are  shown  with  the  stitched  dataset  in  Figure  7.  Investigation  of  Figure  7(d)  at  higher 
magnifications  found  no  discernible  seam  between  the  datasets  and  no  jogs  at  the  grain 
boundaries.  Most  importantly,  the  grains  are  contiguous  across  the  seam  between  the  tiles  and 
are  not  separated  during  the  pixel-grouping/grain  identification  process  permitting  accurate 
measurement  of  grain  morphology  over  large  areas  with  minimal  distortion  if  sufficiently  high 
magnifications  are  used  [Nolze,  2007]. 

In  Figure  8(a),  the  search  regions  from  the  left  and  upper  tiles,  which  were  registered  with  Eqs. 
(15)  -  (18),  are  shown  along  with  the  inverted  L  shaped  search  region  from  the  current  internal 
tile.  The  correlation  between  these  two  arrays  is  shown  in  Figure  8(b)  while  the  Savitzky-Golay 
filter  (linear  fit,  window  size  of  5)  is  shown  in  Figure  8(c).  Background  subtraction,  Figure  8(d), 
essentially  flattens  the  correlation  function  at  this  level  of  resolution,  but  a  single  strong  intensity 
of  -1500  compared  to  the  near  zero  background  marks  the  position  of  maximum  grain  boundary 
overlap.  The  entrywise  product  of  the  search  regions  corresponding  to  this  peak  and  a  higher 
magnification  image  of  the  peak  in  the  cross-correlation  are  shown  in  Figure  9.  The  white  pixels 
remaining  in  the  entrywise  product  reveal  correspond  to  the  grain  boundary  pixels  that  were 
overlapped  at  the  peak  in  the  cross-correlation  function. 

The  complete  dataset  is  shown  in  Figure  10(a)  while  representative  100-pm  x  100-pm  regions 
are  shown  at  higher  magnification  in  Figure  10(b)  -  (d).  These  regions  straddle  locations  where 
four  tiles  meet  and  have  grain  boundaries  with  misorientation  angle  exceeding  5°  highlighted 
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black.  These  images  also  utilize  grayscale  image  quality  overlay  to  highlight  the  seam  between 
the  tiles  in  Figure  10  (b-e).  The  first  two,  b  and  c,  are  typical  examples  of  good  mating  between 
adjacent  tiles  where  there  are  no  obvious  artificial  offsets  created  in  the  grain  boundaries.  The 
arrows  in  Figure  10(d)  call  attention  to  grain  morphologies  that  appear  to  be  incorrect,  including 
a  large  step  in  a  twin  boundary  and  grain  boundary  that  is  precisely  parallel  and  at  the  very  top 
edge  of  the  lower  right  tile.  However,  inspection  of  the  raw  data  reveals  this  was  in  fact  the  true 
morphology  of  the  microstructure  at  this  location.  Figure  10(e),  on  the  other  hand,  shows  the  one 
example  of  imperfect  mating  at  the  seam  that  was  identified  in  the  polycrystalline  dataset.  While 
a  small  ledge  was  created  at  the  seam,  recall  that  this  was  the  position  that  maximized  the 
boundary  overlap  in  the  entire  search  region.  The  effect  of  this  error  on  the  overall  grain 
boundary  and  grain  size  statistics  is  insignificant. 

C.  Two-phase  a  +  P  forged  Ti-6A1-4V 

The  a  +  P  forged  Ti  alloy  dataset,  Figure  4(c),  had  the  finest  grain  size  and  also  used  the  coarsest 
step  size  (Table  1).  Consequently,  there  was  a  higher  density  of  boundaries  in  the  search  region 
and  thus  the  background  level  exceeded  the  peak  associated  with  the  proper  overlap  thereby 
necessitating  background  subtraction.  The  search  regions  for  two  adjacent  tiles,  their  cross¬ 
correlation,  the  estimated  background  (m  =  1,  window  size  =  5)  and  the  background  subtracted 
cross-correlation  are  shown  in  Figure  11.  Following  background  subtraction,  there  is  a  flat 
background  with  an  intensity  of  ~20  while  the  peak  is  just  under  100  giving  a  signal  to  noise 
ratio  of  — 5 : 1 .  Figure  12  shows  the  stitched  tiles  with  an  image  quality  overlay  to  highlight  the 
good  matching  along  the  seams  between  the  tiles.  Despite  the  comparatively  low  signal  to  noise 
ratio,  the  cross-correlation  function  worked  remarkably  well  to  stitch  the  entire  dataset  (Figure 
13).  Due  to  the  large  number  of  grains  in  this  sample,  memory  errors  were  encountered  when 
attempting  to  open  the  complete  stitched  dataset  and  the  final  map  shown  in  Figure  13  was 
created  by  breaking  the  entire  stitched  dataset  into  halves  which  were  able  to  be  opened  in  the 
standard  commercial  packages  to  create  two  crystal  orientation  maps  which  were  merged  in 
standard  image  processing  software. 

V.  Discussion 
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Increased  data  collection  rates  have  necessitated  the  development  of  methods  for  producing  large 
montages  of  data.  In  terms  of  data  collection,  the  generally  accepted  method  now  involves 
collecting  a  large  number  of  individual  tiles  that  each  consist  of  a  standard  beam  raster  scan.  This 
offers  two  possibilities  for  analysis.  The  first  is  to  produce  a  large  amount  of  data  from 
individual  scan  tiles  and  perform  statistical  analyses  of  the  data  in  the  individual  tiles  while  the 
second  is  to  produce  a  single  file  containing  the  combined  infonnation  from  all  of  individual 
scan  tiles.  Each  method  is  useful  in  specific  instances,  but  it  is  generally  desirable  to  have  access 
to  all  of  the  functions  in  commercial  software  for  characterizing  the  state  of  the  microstructure, 
and  so  the  latter  method  was  the  impetus  for  this  work. 

Execution  of  the  source  code  from  its  native  form  is  relatively  fast  on  standard  laptops  and 
desktop  computers  when  Eqs.  (5),  (10)  and  (19)  are  solved  in  the  frequency  domain.  Figure  14 
shows  a  comparison  of  the  processor  time  required  for  various  operations  as  a  function  of  the 
number  of  elements  in  the  search  region  for  a  single  stitching  operation,  i.e.  combining  only  two 
tiles.  The  data  reveal  the  significant  expense  associated  with  performing  correlation  in  the  spatial 
domain  for  any  reasonable  number  of  elements  in  the  search  region.  This  result  supports  the 
choice  of  applying  the  convolution  theorem,  Eq.  (5),  for  performing  automated  stitching.  The 
times  for  this  operation  are  all  less  than  1  second  up  to  -1000  elements.  The  most  costly  function 
in  the  code  is  therefore  identifying  the  pixels  adjacent  to  grain  boundaries,  though  Figure  14 
shows  that  the  operation  is  still  executed  in  a  reasonable  time. 

Prior  to  converting  to  the  frequency  domain  method,  subset  template  matching  was  also 
investigated  as  a  means  for  increasing  speed.  This  method  was  found  to  be  insufficient  in  terms 
of  accuracy  because  the  peak  in  the  correlation  for  a  subset  in  one  location  did  not  always 
correspond  to  the  best  fit  when  the  entire  tile  was  considered.  This  is  because  distortions, 
uncertainties  and  assumptions  made  within  the  SEM  and  EBSD  acquisition  system  [Schwarzer, 
1995;  Nolze,  2007]  lead  to  an  uncertainty  in  the  precise  position  of  grain  boundaries,  particularly 
near  the  periphery  of  scans.  Consider  Figure  7(c)  and  Figure  9  which  show  the  entrywise  product 
of  the  search  regions  at  the  location  that  maximized  the  cross-correlation  function.  The 
boundaries  are  nominally  clear  and  continuous  over  the  entirety  of  the  individual  tiles,  but  this 
clarity  is  only  maintained  in  the  central  region  of  the  entrywise  product  of  these  data  (Figure 
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7(c)).  The  primary  cause  of  this  effect  is  the  trapezium  (trapezoidal)  distortion  that  is  prevalent  at 
low  magnifications  where  large  beam  deflection  angles  must  be  used  to  reach  the  periphery  of 
the  scan  [Nolze,  2007],  In  practice,  one  usually  adjusts  the  focus  to  be  best  at  the  position  Ysizel2, 
which  also  corresponds  to  the  location  where  the  mismatch  between  two  adjacent  tiles  caused  by 
trapezium  distortion  is  the  smallest.  Thus,  it  is  tempting  to  use  this  as  a  template  centered  on 
YsizJ2  to  perfonn  the  cross-correlation.  However,  it  is  not  trivial  to  ensure  the  sample  surface  is 
parallel  with  the  x-y  plane  of  the  SEM  stage  and  the  act  of  moving  the  sample  may  cause  the 
distance  between  the  pole  piece  and  the  sample  surface  to  change  which  changes  the  y- 
coordinate  in  the  scan  at  which  the  best  focus  is  achieved.  Said  another  way,  if  the  sample 
surface  is  not  truly  at  the  70°  tilt  angle  that  is  expected  by  the  EBSD  acquisition  system,  this 
introduces  yet  another  distortion  [Nolze,  2007].  This  distortion,  which  either  foreshortens  or 
lengthens  the  features  of  the  microstructure  depending  on  the  deviation  from  70°,  will  not 
necessarily  result  in  obvious  errors  at  tile  seams  as  it  will  be  similar  for  adjacent  tiles  in  a  given 
dataset.  However,  it  is  important  to  recognize  that  this  distortion  exists  as  it  may  skew 
measurements  of  specific  features  of  the  microstructure. 

There  is  also  reduced  beam  positioning  sensitivity  at  large  deflections  from  the  optic  axis  which 
adds  uncertainty  to  the  precise  location  of  the  beam  when  a  particular  data  point  was  collected. 
Also,  because  the  beam  is  slightly  defocused  at  the  longer  or  shorter  working  distance 
corresponding  to  the  top  and  bottom  of  the  scan  tile,  respectively,  grain  boundaries  usually 
appear  more  diffuse  in  image  quality  or  band  contrast  maps  at  these  locations.  Since  the  present 
grain  boundary  fonnulation  operates  on  a  pixel-to-pixel  basis,  it  is  easy  to  see  how  the 
aforementioned  problems  may  cause  some  parts  of  grain  boundaries,  or  entire  boundaries 
altogether,  to  not  appear  in  the  entry- wise  product  of  the  grain  boundary  arrays  corresponding  to 
the  maximum  intensity  in  the  cross-correlation  (Figure  6(c)  and  Figure  7(c)  and  Figure  9).  In 
practice,  the  effect  of  SEM-induced  distortions  can  be  minimized  by  simply  using  higher 
magnifications.  Nolze  [Nolze,  2007]  has  shown  that  a  distortion  of  approximately  2%  exists  at 
magnifications  of  400x  while  the  1%  threshold  is  crossed  at  a  magnification  of  700x  for  the 
instrument  used  in  the  investigation.  In  our  SEM  (FEI XL30  ESEM),  500x  magnification  at  25 
mm  working  distance  while  employing  dynamic  focus  routines  was  sufficient  for  multi-tile 
mapping  with  minimal  distortion  observed  between  the  tilt  corrected  EBSD  data  and  images 
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taken  at  zero  degree  tilt.  Imaging  of  transmission  electron  microscopy  grids  at  zero  tilt  and 
inclined  tilt  is  recommended  for  quantifying  image  distortions  inherent  to  each  users  SEM. 

In  cross-correlation  of  images  with  a  distribution  of  grayscale  intensities,  it  is  common  practice 
to  scale  the  intensity  values  to  an  interval  centered  on  zero  so  that  overlapping  local  minima  in 
an  image  pair  contribute  toward  identifying  a  correct  match,  in  addition  to  overlapping  local 
maxima.  In  the  present  formulation,  this  is  analogous  to  setting  grain  interiors  to  -1,  while  still 
setting  boundaries  to  +1,  prior  to  performing  cross-correlation  via  Eqs.  (5),  (10)  and  (19).  In 
preliminary  tests,  this  was  found  to  decrease  the  strength  of  the  peak  signal  relative  to  the 
background  level.  More  importantly,  however,  the  strength  of  the  peak  in  the  cross-correlation 
decreased  relative  to  its  immediate  surroundings  and  this  trend  was  consistent  for  most  datasets. 
Ultimately,  the  same  solution  was  reached  after  background  subtraction.  Using  the  [-1,1]  range 
merely  meant  that  every  dataset  needed  background  subtraction  leading  to  overall  longer 
execution  times  on  datasets  with  a  low  density  of  boundaries.  Thus,  the  final  version  of  software 
uses  data  on  the  interval  [0,1]  and  the  use  of  background  subtraction  is  optional.  If  enabled,  the 
user  is  required  to  enter  the  order  of  the  polynomial  and  the  window  size.  In  most  cases,  a  simple 
linear  fit,  m  =  1  in  Eq.  (24),  over  3-7  data  points  obtained  a  sufficient  background. 

In  closing,  we  recognize  that  in  order  to  function  properly,  the  present  formulation  requires 
EBSD  data  of  relatively  high  quality.  This  implies  that  there  are  very  few  improperly  indexed 
pixels.  A  few  guidelines  for  data  collection  were  developed  while  writing  and  testing  the  code: 

If  possible,  prepare  the  sample  by  electropolishing.  This  is  among  the  most  forgiving 
ways  to  polish  metallic  materials  for  large-area  EBSD  evaluation. 

The  sample  surface  should  be  carefully  aligned  with  the  X-Y  translation  plane  of  the 
stage  either  by  precise  parallel  polishing  of  the  sample,  or  through  the  use  of  an 
adjustable  sample  holder. 

The  SEM  magnification  should  be  sufficiently  high  that  obfuscation  of  grain  boundaries 
by  SEM-induced  spatial  distortions  are  minimized. 

The  dynamic  focus  feature  should  be  enabled,  if  available. 

The  size  of  the  overlap  region  assigned  during  data  collection  is  a  function  of  the  step 
size,  grain  size  and  grain  morphology.  Users  will  develop  a  feel  for  this  for  their  own 
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samples  to  obtain  high  quality  cross-correlations.  [Humphreys,  2001],  however,  has 
suggested  that  a  minimum  of  8-pixels  per  constituent  is  necessary  to  obtain  appropriate 
size  statistics  and,  as  a  minimum.  In  the  present  case,  a  minimum  of  about  8  pixels  per 
grain  appears  to  be  sufficient  to  apply  these  methods. 

V.  Summary 

In  summary,  we  have  demonstrated  a  suitable  method  for  automatically  combining  multi-tile 
EBSD  datasets  with  high  precision  to  maintain  precise  grain  morphology  in  single  and  two-phase 
polycrystalline  materials.  The  algorithm  simply  maximizes  the  amount  of  grain,  phase  and  (or) 
grain  and  phase  boundary  overlap  between  adjacent  tiles  in  EBSD  scans  and  therefore  should  be 
applicable  to  more  complex  materials  as  well.  Equations  were  presented  for  each  of  the  three 
possible  tile  positions:  first  row  of  data,  first  column  of  data  and  internal  tiles.  Standard  cross¬ 
correlation  algorithms  were  successfully  applied  to  recover  the  linear  shifts  necessary  to  translate 
the  EBSD  data  into  proper  registration  with  no  artificial  seams  between  the  identified 
boundaries.  In  datasets  with  dense  boundaries,  background  estimation  was  performed  using  a 
Savitzky-Golay  moving  window  filter  row-wise,  then  column  wise,  and  this  result  was 
subtracted  from  the  cross-correlation.  In  most  cases,  a  simple  linear  filter  in  two  directions  was 
sufficient. 
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(a)  (b)  (c)  (d)  (e)  (f) 

Figure  1 .  (a)  and  (b)  Binary  representation  of  grain  boundaries  in  the  overlapping  region  of  two 
adjacent  tiles  of  a  polycrystalline  nickel-base  superalloy  EBSD  dataset,  (c)  The  cross-correlation 
of  (a)  and  (b).  (d)  Horizontally  and  vertically  mirrored  version  of  (b).  (e)  The  cross-correlation  of 
(a)  and  (d),  which  is  equivalent  to  the  convolution  of  (a)  and  (b).  (f)  The  instantaneous  area  of  the 
overlapping  regions  of  the  tiles  at  each  pixel  in  the  cross-correlation  and  convolution.  The  insets 
in  (c)  and  (e)  show  the  regions  around  the  peak  value  in  each  image  at  a  magnification  of  300 
times. 
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Figure  2.  Schematic  diagram  of  the  various  cases  encountered  during  stitching,  (a)  First  row,  (b) 
first  tile  of  a  new  row  of  data,  (c)  internal  tiles  and  the  associated  spatial  (x,y)  and  EBSD  (X,Y) 
reference  frames.  In  (d),  the  left  and  upper  tiles  for  the  internal  tile  case  have  been  stitched  using 
the  previously  calculated  values  prior  to  performing  cross-correlation  with  the  search  region  of 
the  current  tile. 
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Figure  3.  Grain  boundary  detection  using  [2x2]  kernel.  The  individual  components  are  shown  in 
(a),  (b)  and  (c)  while  the  combined  set  appears  in  (d). 
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Figure  4.  Crystal  orientation  maps  for  the  first  (left)  and  second  (current)  tiles  of  (a)  coarse 
grained  a-Titanium,  (b)  a  2x2  array  of  a  hot  pack  rolled  and  recrystallized  Ni-Al-Cr  superalloy 
and  (c)  a  +  (3  processed  Ti-6A1-4V.  The  scale  bars  are  (a)  20  pm,  (b)  30  pm  and  (c)  50  pm  and 
show  the  range  covered  by  the  search  region  in  each  tile. 
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Figure  5.  Grain  boundaries  in  the  (a)  left  and  (b)  current  tiles,  (c)  Cross-correlation  of  (a)  and  (b). 
(d)  Savitzky-Golay  smoothing  of  (c)  performed  with  a  first  order  (linear)  polynomial  and  a 
window  size  of  3.  (e)  Difference  between  (c)  and  (d). 
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(a)  (b)  (c) 


Figure  6.  Single  phase  Ti-7A1  (wt%).  Binary  data  showing  the  grain  boundaries  in  only  the 
overlapping  parts  of  the  (a)  left  and  (b)  current  tiles,  and  the  (c)  entry-wise  product  of  (a)  and  (b) 
when  the  cross-correlation  function  is  maximized,  (d)  Crystal  orientation  map  colored  according 
to  the  standard  triangle  for  hexagonal  materials  shown  in  the  lower  right  comer  of  the  image. 
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Figure  7.  y+  y’  strengthened  Ni-Al-Cr-base  superalloy.  Binary  representation  of  the  grain 
boundaries  in  only  the  overlapping  regions  of  the  (a)  left  tile  and  (b)  current  tile,  and  the  (c) 
entry-wise  product  of  (a)  and  (b)  when  the  cross-correlation  function  is  maximized,  (d)  Crystal 
orientation  map  colored  according  to  the  standard  triangle  for  cubic  materials  shown  in  the  lower 
right  corner  of  the  image.  These  two  tiles  were  combined  using  the  correlation  shown  in  Figure 
1(c). 
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Figure  8.  (a)  Grain  boundaries  in  the  search  region  of  the  combined  left  and  upper  tiles  and  the 
current  internal  tile  in  the  polycrystalline  nickel  sample,  (b)  Correlation  of  the  two  regions  shown 
in  (a),  (c)  Savitzky-Golay  estimation  of  the  background,  linear  fit  over  a  window  size  of  5.  All 
data  are  shown  in  pixels  where  1  pixel  =  0.5  pm.  (d)  Background  corrected  cross-correlation. 

The  peak  intensity  is  shown  at  higher  magnification  in  Figure  9. 
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Figure  9.  Entry-wise  product  of  the  arrays  in  Figure  8(a)  when  the  background  corrected 
correlation,  Figure  8(d),  is  maximized.  The  boxed  region  around  the  peak  in  Figure  8(d)  is  shown 
at  higher  magnification. 
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Figure  10.  (a)  Normal  direction  inverse  pole  figure  map  with  image  quality  overlay  and  grain 
boundaries  >5°  colored  black  for  the  42-tile,  polycrystalline  Ni  superalloy  data  set.  (b),  (c)  and 
(d)  are  100-pm  x  100-pm  regions  within  (a)  shown  at  a  magnification  of  400x.  Faint  vertical 
bands  in  image  quality  highlight  the  seams  between  adjacent  tiles.  In  (d),  some  flat  segments  of 
grain  boundaries  parallel  to  the  seam,  or  large  ledges  in  twins  give  the  impression  of  an  error,  but 
are  actually  present  in  the  raw  data.  In  (e),  the  arrows  show  an  example  of  an  artifical  step 
created  by  the  stitching  process. 
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Figure  11.  Grain  and  phase  boundaries  in  the  search  regions  of  the  (a)  left  and  (b)  current  tiles  of 
an  a  +  P  forged  Ti-6A1-4V  alloy,  (c)  correlation  of  (a)  and  (b).  (d)  Linear  Savitzky-Golay  filter 
(m  =1)  with  a  window  size  of  5.  (e)  The  difference  between  (c)  and  (d). 
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Figure  12.  Crystal  orientation  map  (see  legend  in  Figure  6)  of  a  +  P  forged  Ti-6A1-4V  alloy 
combined  with  the  correlation  shown  in  Figure  11.  The  image  quality  overlay  is  used  to  highlight 
the  seam  between  the  tiles  which  is  in  approximately  the  center  of  the  image. 
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Figure  13.  Crystal  orientation  map  (see  legend  in  Figure  6)  of  the  entire  dataset  for  the  a  +  p 
forged  Ti-6A1-4V  alloy. 
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Figure  14.  Processor  time  as  a  function  of  number  of  elements  in  the  search  region  for  the  most 
computationally  expensive  calculations  in  the  code. 
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